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This  report  investigates  the  plane  stress  capability  of  the  FEARS  (Finite 
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the  results,  and  provide  the  conclusion. 
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ABSTRACT 


This  report  investigates  the  plane  stress 
capability  of  the  FEARS  (Finite  Element  Adaptive 
Refinement  Solver)  program  which  features  both 
adaptive  meshing  and  a  posteriori  error  bounding. 

Two  benchmark  plane  stress  problems  were  solved  by 
both  the  FEARS  and  NASTRAN  (NASA  Structural  Analysis) 
programs  for  comparison.  This  report  presents  only 
the  basic  computational  results.  Other  reports  to 
follow  shortly  will  motivate  the  experimental  proce¬ 
dures,  analyze  the  results,  and  provide  the  conclu¬ 
sion. 


ADMINISTRATIVE  INFORMATION 

This  work  was  performed  under  the  DTNSRDC ' s  Independent  Research  Program, 
Program  Element  6U52N,  Task  Area  ZR0140201,  DTNSRDC  Work  Unit  1844-140.  A 
contract  arrangement  was  entered  into  with  Prof.  Babuska  of  the  Institute  for 
Physical  Science  and  Technology,  University  of  Maryland. 

1.  INTRODUCTION 

This  report  presents  the  basic  results  of  the  numerical  solutions  of  two 
benchmark  problems  in  plane  elasticity.  The  purpose  of  this  computation  was 
to  compare  the  performance  of  the  NASTRAN  (NASA  Structural  Analysis)  and  FEARS 
(Finite  Element  Adaptive  Refinement  Solver)  programs.  A  subsequent  report 
will  analyze  the  results  and  present  conclusions. 

Section  2  of  this  report  describes  the  two  benchmark  problems.  Section  3 
briefly  describes  the  use  of  the  NASTRAN  program  to  solve  the  two  problems, 
and  also  describes  the  data  generator  program,  DGNEW,  used  to  prepare  data 
cards  for  NASTRAN.  Section  4  presents  the  results  obtained  by  NASTRAN  for 
both  benchmark  problems.  Section  5  briefly  describes  the  computation  with  the 
FEARS  program  and  Section  6  presents  the  basic  results  of  the  computation. 


2.  THE  BENCHMARK  PROBLEMS 

The  two  plane  stress  benchmark  problems  are  illustrated  in  Figures  1  and 
2.  Figure  1  shows  the  domain  and  the  load  on  the  sides  of  Problem  1.  Because 
the  domain  is  symmetric  only  the  upper  right  hand  shaded  portion  need  be 
considered.  We  shall  give  the  values  of  the  parameters  later. 
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Problem  2  uses  the  same  configuration  as  problem  1  but  with  a  crack  In 
the  vertical  direction.  The  length  of  the  crack, X,  is  -|<h-r).  See  Figure  2. 

A  detailed  justification  for  the  selection  of  these  two  problems  will  be  given 
in  a  forthcoming  report.  The  present  investigation  discusses  the  FEARS 
program  with  respect  to  the  presence  of  singularities  in  the  geometry  and  the 
solution.  Problem  1  clearly  degenerates  as  r  approaches  h,  and  Problem  2  has 
a  singular  solution  at  the  tip  of  the  crack.  Both  Problems  1  and  2  are  plane 
stress  problems  where  E,  Young's  modulus  of  elasticity,  is  3.0  x  10^  and  v, 
Poisson's  ratio,  Is  0.3.  These  two  problems  are  discussed  individually  in 
more  detail . 


Problem  1. 

The  following  parameters  were  used 
h  ■  1 
l  •  6 

p  ■  1  (normal  stress) 

The  solution  to  problem  1  consists  of  the  displacement  vector 


(u,v)T 

and  the  tensor  of  stresses 


T 

The  problem  is  to  find  the  solution  (u,v)  with  sufficient  accuracy  with 

respect  to  the  energy  norm  and  the  stress  concentration  factor  Z  of  the  stress 

a  where 
xx 


Z 


max  o 
_ xx 

£h_ 


h-r 


Problem  2. 

Problem  2  has  the  same  geometry  as  Problem  1,  but  with  a  crack  of  length 
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X.  The  following  parameters  were  used 

h  -  1 
1-6 

p  -  1  (normal  stress) 
r  -  .7 
X  -  .15 

T 

The  problem  Is  to  find  the  solution  (u,v)  and  the  stress  Intensity 
factor,  Kj ,  with  sufficient  accuracy.  is  the  coefficient  of  the  principal 
singular  part  of  the  solution. 

As  before  the  problem  involves  the  energy  e  and  a  constant,  K  ,  the 
stress  Intensity  factor,  defined  to  be  the  coefficient  of  the  first  singular 
term  of  the  solution  at  the  tip  of  the  crack. 

For  more  information  on  the  following  formulas  the  reader  is  referred  to 
Pu ,  Hussain ,  et  al .  [ 2 ]  and  chapter  2  of  Morosov  and  Nikischov  [ 3 ] .  The 
solution  is  singular  at  the  tip  of  the  crack.  If  we  Introduce  polar  coordi¬ 
nates  as  shown  in  Figure  3  we  obtain 


Figure  3 .  Polar  €e irdlnate  Scheme  for  Problem  2 
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u 


(1) 


F  (0) 
u 


v  -  — W—  F  (0) 

2G12II  V  ' 

where 

Fu<0)  -  «in|(ie  +  1  -  2cos2|) 

Fv(Q)  -  cosj(ic  -  1  +  2sin2  y) 
E 

and  G  ■  +  v)  is  the  shear  modulus. 

In  the  case  of  plane  stress 


K 


3  -  y 
1  +  V 


(2) 

(3) 

(4) 


These  formulas  allow  us  to  compute  In  different  ways.  If  0  is  (0,0) 
and  A  ■  (Lsin  0,Lcos  0)  then  we  can  appropriate  by 


or 


(5a) 


(2GV£nrv(A)  -  v(0)1 


VTfv(©) 


Here  (u(0),v(0))  is  the  displacement  of  the  tip  of  the  crack  and  (u(L),v(L)) 
is  the  displacement  of  the  grid  point  A  located  a  distance  L  from  the  tip  of 
the  crack  and  where  the  displacement  vector  makes  an  angle  0  with  the  y-axis, 
as  shown  in  Figure  3. 
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Alternatively,  another  term  can  be  added  to  the  expansion  of  (u,v),  that 
is 


+  C2VT  +  C3r 

(6a) 

+  C,  VT  +  C  r 

J  *  0 

(6b) 

where  the  are  functions  of  0.  If  two  collinear  grid  points  make  an  angle  0 
with  the  y-axis,  and  they  are  located  at  distances  L  and  £L  (where  0<£<1) 


(7a) 


or 


x  -(2G^H)Cv  (£L)  -  Sv(L)  -  (1-0  v(0) ] 

1  (VT-\)>/rFv(0) 

(7b) 


Finally  can  also  be  computed  by  the  energy  method.  When  the  energy 
C  of  the  solution  is  a  function  of  the  length  of  the  crack  X  then 


(8) 


However,  this  computation  requires  a  second  finite  element  model  where  the 
crack  length  has  been  changed  by  an  amount  dX. 


3.  DESCRIPTION  OF  NASTRAN  COMPUTATION 
A  plane  stress  problem  is  input  to  NASTRAN  via  the  ’’bulk  data  deck" 
consisting  of  the  grid,  connection,  force,  constraint,  and  material  properties 
cards.  The  geometry  of  the  problem  is  specified  on  the  grid  and  connection 
cards  which  must  adhere  to  a  rigid  format  decreed  by  NASTRAN.  Since  preparing 
these  cards  by  hand  is  both  tedious  and  laborious  a  computer  program,  DGNEW, 
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(a  data  generator)  was  written  to  generate  the  mesh  from  which  these  grid  and 
connection  cards  are  produced.  Often  the  data  generator  is  a  fairly  complicated 
program  in  its  own  right  as  in  the  present  case. 

The  data  generator  DGNEW  generates  the  meshes  for  both  problems  1  and  2. 

The  parameters  N^,  ,  r,  h,  £,  a,  and  KK  are  read  from  tape  5  in  free  format. 

Nj  is  the  number  of  partitions  of  side  AB  of  subdomain  (Figure  Ab) ,  and  Nj 
is  the  number  of  partitions  of  side  HG  of  subdomain  0^  (Figure  Ad),  r  is  the 
radius  of  the  arc  of  the  circle  in  subdomain  0^  (Figure  Aa) .  h  is  the  length 
of  side  HG  of  subdomain  (Figure  Aa) .  I  is  the  length  of  side  BDH  (Figure 
Aa).  a  is  the  desired  aspect  ratio.  If  the  value  of  the  last  parameter  KK  is 
zero  then  it  is  assumed  there  is  no  crack.  Otherwise  the  crack  begins  at  the 
node  specified  by  the  value  of  KK.  DGNEW  generates  card  images  of  the  appro¬ 
priate  grid,  connection,  force,  and  constraint  cards.  A  listing  of  DGNEW 
which  utilizes  quarter  points  for  the  crack  tip  of  problem  2,  is  given  in  the 
appendix  to  this  report . 

The  NASTRAN  data  deck  also  contains  control  cards  pertaining  to  the 
computation  and  to  the  printing  of  output.  These  cards  must  also  adhere  to  a 
specified  NASTRAN  format. 

The  NASTRAN  IS2D8  element  was  used  for  these  computations .  It  is  the 
usual  two-dimensional,  quadratic,  isoparametric,  plane  stress  element  with 
eight  nodes,  a  so-called  "serendipity"  element.  Stresses  at  the  nodes  are 
extrapolated  from  stress  computed  on  3  x  3  array  of  Gauss  integration  points . 
Elements  of  this  type  are  of  degree  2,  that  is,  for  a  smooth  solution  and 
uniform  meshes  whose  elements  are  squares  with  side  h  the  rate  of  convergence 
in  the  energy  norm  is  O(h^)  *  0^—^ where  N  is  the  number  of  degrees  of  freedom. 

For  Problem  1  only  elements1  of  this  "serendipity"  type  were  used.  For 
Problem  2  two  different  meshes  were  used.  The  first  mesh  consists  entirely  of 
"serendipity"  elements  as  described  above.  The  second  mesh  consists  of  the 
same  "serendipity"  elements  except  for  those  two  elemets  which  have  had  the 
"midpoints"  of  their  two  sides  adjacent  to  the  vertex  (which  is  the  tip  of  the 
crack)  changed  into  "quarter  points".  That  is,  on  the  two  sides  of  these  two 
elements  which  intersect  at  the  tip  of  the  crack,  the  "midpoints"  are  now 
located  only  a  quarter  of  a  length  of  their  side  away  from  the  tip  of  the 
crack,  instead  of  half  the  length  of  their  side.  We  shall  refer  to  these 
elements  as  "quarter-point  'serendipity’"  elements. 
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The  mesh  generator  constructs  a  mesh  for  ft  using  the  physical  domensions 
(r,jt,h)  of  the  domain  ft  pictured  in  Figure  4a,  the  prescribed  number  of  left 
and  right  partitions,  (Figure  4b)  and  (Figures  4c  and  4d),  and  the 
specified  aspect  ratio  a  for  the  elements.  In  the  course  of  generating  this 
mesh,  two  other  integers  (N^,  N^)  and  a  real  number  q  (r<q<&)  are  computed 
from  this  input  data.  is  the  number  of  partitions  of  the  side  BD  of 

subdomain  r^  (Figure  4b).  is  the  number  of  columns  of  elements  In  sub- 
domain  R^  (Figure  4d).  q  is  the  magnitude  of  side  BD  of  subdomain  ftj.  These 
nine  parameters  (r,  £,  h,  ,  N^,  a,  N^,  N^.  and  q)  completely  describe  the 
mesh  constructed  for  ft.  The  original  domain  ft  is  then  subdivided  into  three 
subdomains  R^,  R^,  and  R^  as  shown  in  Figure  4a.  The  subdomains  are  ft^  * 
ABCD,  R2  -  CDFE,  and  ft3  -  FBHG. 

The  mesh  of  R^  is  shown  in  Figure  4b.  The  side  BD  is  divided  into 

intervals,  thus  determining  N^+l  grid  points.  (The  method  for  determining 

and  q  is  described  later  in  this  section.)  The  arc  AC  is  divided  into 

non-uniform  Intervals  determined  by  the  N,  angles  1 1». ,  i and  has  N.+l 

J  1  /  N3  j 

grid  points  along  it  determined  by  the  intervals.  The  respective  grid  points 
on  BD  and  AC  are  connected  by  straight  lines  which  in  turn,  are  divided  into 
Nj  equal  Intervals  thus  determining  N^+l  grid  points  along  each  line.  The 
mesh  of  ft^  is  determined  by  these  (Nj+1)(N3+1)  grid  points. 

The  angles  and  are  chosen  such  that  the  aspect  ratio  of 

any  element  is  approximately  a.  The  procedure  is  as  follows: 


where  i is  the  approximate  length  of  the  line  segment  connecting  the 
respective  grid  points  on  AC  and  BD.  is  given  by  the  formula 
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Figure  4a  -  Partition  of  Domain  fi 


q  N3  elements 


Figure  4b  -  The  Mesh  of  Subdomain  0. 


Figure  4c  -  The  Mesh  of  Subdomain 


Figure  4d  -  The  Mesh  of  Subdomain 
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where  i  -  CD  is  defined  by 


(r  +  &)cos$  ■  h 
and 

!•  ■  AB  "  h-r. 

o 

If  T.  Is  defined  as  ij/.  +  +...$'  then  for  some  integer  j  there  will 

w  1  >  1  *  1  ° 

which  will  be  — (3  where  $  is  the  radial  angle  of  subdomain  £1^  In  Figure  4b. 
If  is  the  smallest  Integer  for  which  this  is  true,  the  angles  <1^  are 
obtained  by  scaling  down  the  ^  thus 

*i  - 

The  mesh  in  is  shown  in  Figure  4c.  The  sides  EC  and  FD  are  divided 
into  equal  segments .  The  respective  grid  points  on  each  side  are  connected 
by  straight  lines  each  of  which,  in  turn,  is  divided  into  equal  parts, 
giving  (Nj+lXN^+l)  grid  points  for  the  mesh  of  The  required  aspect  ratio 

of  the  elements  in  ^determines  q.  The  angle  Y  is  obtained  by  solving 

tan  y  - 

The  mesh  of  is  shown  in  Figure  4d.  N^,  the  number  of  equal  partitions 

of  both  sides  FG  and  DH,  is  determined  by  the  aspect  ratio  a  thus 


A  -  ■££=& 

*4 

Since  the  sides  FD  and  GH  are  both  divided  into  equal  parts ,  there  are 
(N^+1)(N2+1)  grid  points  for  03  which,  in  turn,  determine  a  uniform 
rectangular  mesh  for 

The  elements  determined  by  the  grid  points  of  the  mesh  of  £1  are 
isoparametric  with  eight  grid  points;  the  mid-points  are  located  in  the  middle 
of  the  sides.  When  the  side  of  an  element  is  a  circular  arc,  arc  length  is 
used  to  locate  the  mid-point. 
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Table  1  presents  the  various  aeshes  (number  of  elements  and  grid  points) 
generated  for  one  value  of  r  (.999)  with  different  aspect  ratios  and  several 
choices  of  Nj  and  to  show  the  dependence  of  the  number  of  elements  and  DOF 
on  the  parameters.  The  Table  1  colu^a  headings  are  identified  as  follows: 

1  is  the  line  number  of  the  table 

K  is  the  microfiche  identification  number  for  the  particular  NASTRAN  run 
a  is  the  aspect  ratio 

Nj  is  the  number  of  the  left  hand  partitions 
N2  Is  the  number  of  the  right  hand  partitions 

is  the  number  of  elements  in  the  first  subdomain 

1^2  is  the  number  of  elements  in  the  second  subdomain 

Nq ..  is  the  number  of  elements  in  the  third  subdomain 

Nft  is  the  total  number  of  elements  in  the  domain 
GP  is  the  total  number  of  grid  points  (nodes) 

DOF  is  the  number  of  degrees  of  freedom 

TABLE  1  -  MESHES  GENERATED  FOR  R  -  0.999  IN  PROBLEM  1 


I 

K 

a 

N1 

N2 

Nfil 

n«2 

% 

GP 

DOF 

1 

8 

1.5 

4 

3 

96 

12 

24 

132 

475 

916 

2 

9 

1.5 

6 

4 

198 

24 

44 

266 

907 

1766 

3 

10 

1.5 

8 

5 

336 

40 

70 

446 

1477 

2892 

4 

17 

1.0 

12 

8 

1068 

96 

264 

1428 

4569 

9022 

5 

18 

1.0 

8 

5 

480 

40 

100 

620 

2047 

4020 

6 

19 

1.0 

4 

3 

132 

12 

362 

186 

645 

1248 

7 

20 

.75 

12 

8 

1188 

96 

352 

1636 

5235 

10332 

8 

21 

.5 

8 

5 

792 

40 

200 

1032 

3401 

6688 

9 

22 

.3 

8 

5 

792 

40 

330 

1162 

3843 

7520 

Figure  5  shows  a  sample  mesh  for  Problem  1. 


Figure  5  -  A  Simple  Mesh  for  Problem  1 
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4.  NASTRAN  RESULTS 

In  the  solution  of  the  two  benchmark  problems  using  NASTRAN,  as  shown 
schematically  in  Figures  1  and  2,  h  ■  1,  i  ■  6,  p  ■  1,  and  r  was  varied.  The 
results  for  Problem  1  are  given  In  Table  2.  The  meanings  of  the  column 
headings  are  as  previously  defined  with  the  following  additions: 

R  is  the  radius  of  the  circle  triiose  arc  is  a  side  of  domain  ft. 

SCF  is  the  stress  concentration  factor  £  previously  defined  in  the 

text. 

CP  is  the  CYBER  7400  central  processor  time.  (The  CYBER  7400  is 
essentially  a  CDC  6600.) 

The  asterisks  on  some  of  the  stress  concentration  factors  in  Table  2  mean 
that  the  maximum  stress  did  not  occur  at  the  expected  place  (0,r)  for  that 
NASTRAN  solution. 

The  scheme  of  Problem  2  is  given  in  Figure  2,  and  the  parameters  used 
were  h*  1,1*6,  p  ■  1,  r  ■  .7  and  A  -  .15.  As  noted  before.  Problem  2  was 
solved  in  two  different  ways,  that  is,  with  and  without  quarter-point 
elements.  These  results  are  presented  in  Tables  3a  (all  "serendipity" 
elements  -  no  quarter-point  elements)  and  3b  (all  "serendipity"  elements 
except  for  the  two  elements  that  have  the  tip  of  the  crack  for  a  vertex  - 
these  are  quarter-point  elements).  The  column  heading  KK  refers  to  the  grid 
point  (node)  at  which  the  tip  of  the  crack  is  located.  i8  t*ie  maximum 

stress  which  is  located  at  the  tip  of  the  crack.  The  other  headings  were 
previously  defined.  It  should  be  noted  that  the  exact  stresses  are  infinite 
at  the  tip  of  the  crack.  Table  3c  presents  the  results  when  the  tip  of  the 
crack  for  run  3S  in  Table  3b  is  perturbed  from  (0.,.85)  to  (0. , .868750) . 

As  stated  in  Section  2,  the  stress  intensity  factor  Kj  can  also  be 
computed  by  the  energy  release  method  of  Equation  8.  To  make  use  of  this 
procedure  the  y-coordinate  of  the  tip  of  the  crack  was  perturbed  from  y  *  .85 

(run  3S)  to  .868750  (run  4S)  and  Problem  2  was  solved  with  NASTRAN  for  N^  * 

16,  Nj  ■  12.  This  perturbation  of  the  y-coordinate  was  easily  effected  by 
taking  the  value  of  the  KK  parameter  of  the  data  generator  DGNEW  to  be  19, 
thus  putting  the  tip  of  the  crack  on  the  extreme  left  hand  grid-point  common 
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TABLE  2  -  PROBLEM  1  RESULTS 


I 

K 

R 

N1 

N2 

a 

SCF 

SE  X  10 

GP 

DOF 

CP 

1 

61 

0.05 

8 

6 

1.0 

2.431 

1.0009 

899 

266 

1694 

206 

2 

62 

0.05 

16 

12 

1.0 

2.743 

1.0009 

3343 

1048 

6480 

974 

3 

55 

0.1 

4 

3 

1.0 

2.312 

1.0037 

252 

67 

452 

70 

4 

56 

0.1 

8 

6 

1.0 

2.624 

1.0039 

879 

260 

1656 

202 

5 

57 

0.1 

16 

12 

1.0 

2.746 

1.0040 

3267 

1024 

6332 

951 

6 

52 

0.2 

4 

3 

1.0 

2.437 

1.0162 

241 

64 

6432 

68 

7 

53 

0.2 

8 

6 

1.0 

2.545 

1.0164 

859 

254 

1618 

197 

8 

54 

0.2 

16 

12 

1.0 

2.544 

1.0165 

3241 

1016 

6284 

941 

9 

49 

0.3 

4 

3 

1.0 

2.384 

1.0391 

244 

65 

440 

68 

10 

50 

0.3 

8 

6 

1.0 

2.391 

1.0394 

845 

250 

1594 

195 

11 

51 

0.3 

16 

12 

1.0 

2.373 

1.0394 

3227 

1012 

6262 

934 

12 

46 

0.4 

4 

3 

1.0 

2.286 

1.0761 

244 

65 

440 

68 

13 

47 

0.4 

8 

6 

1.0 

2.264 

1.0764 

877 

260 

1660 

201 

14 

48 

0.4 

16 

12 

1.0 

2.254 

1.0764 

3301 

1036 

6414 

959 

15 

37 

0.5 

4 

3 

1.0 

2.198 

1.1339 

258 

69 

468 

70 

16 

38 

0.5 

8 

6 

1.0 

2.184 

1.1341 

903 

268 

1712 

206 

17 

39 

0.5 

16 

12 

1.0 

2.177 

1.1341 

3451 

1084 

6714 

1013 

18 

40 

0.6 

4 

3 

1.0 

2.145 

1.2247 

261 

70 

476 

71 

19 

41 

0.6 

8 

6 

1.0 

2.134 

1.2249 

935 

278 

1778 

214 

20 

42 

0.6 

16 

12 

1.0 

2.129 

1.2249 

3575 

1124 

6966 

1009 

21 

32 

0.7 

4 

3 

1.0 

2.107 

1.3746 

289 

78 

532 

76 

22 

44 

0.7 

8 

6 

1.0 

2.100 

1.3748 

993 

296 

1896 

227 

23 

45 

0.7 

16 

12 

1.0 

2.097 

1.3748 

3749 

1180 

7318 

1124 

24 

1 

0.8 

4 

3 

1.5 

2.078 

1.6499 

234 

63 

432 

72 

25 

23 

0.8 

4 

3 

1.0 

2.076 

1.6501 

306 

83 

568 

88 

26 

2 

0.8 

6 

4 

1.5 

2.074 

1.6503 

421 

120 

792 

126 

27 

24 

0.8 

8 

6 

1.0 

2.072 

1.6504 

1077 

320 

2066 

307 

28 

25 

0.8 

16 

12 

1.0 

2.071 

1.6505 

4073 

1284 

7970 

1796 

29 

26 

0.85 

4 

3 

1.0 

2.061 

1.8943 

320 

87 

596 

91 

30 

27 

0.85 

8 

6 

1.0 

2.058 

1.8948 

1155 

346 

2222 

327 

31 

28 

0.85 

16 

12 

1.0 

2.057 

1.8948 

4335 

1368 

8496 

1919 

32 

3 

0.9 

4 

3 

1.0 

2.045 

2.3146 

348 

95 

652 

101 

33 

14 

0.9 

6 

4 

1.5 

2.044 

2.3149 

467 

134 

886 

139 

34 

29 

0.9 

8 

6 

1.0 

2.043 

2.3155 

1239 

372 

2392 

353 

35 

30 

0.9 

16 

12 

1.0 

2.042 

2.3156 

4647 

1468 

9122 

2105 

1A 


TABLE  2  (continued) 


I 

K 

R 

N1 

*2 

a 

SCF 

SE  x  10' 

GP 

DOF 

CP 

36 

4 

0.95 

4 

3 

1.5 

2.024 

3.2750 

293 

80 

552 

86 

37 

31 

0.95 

4 

3 

1.0 

2.024 

3.2794 

404 

111 

764 

110 

38 

13 

0.95 

6 

4 

1.5 

2.024 

3.2802 

547 

158 

1046 

162 

39 

32 

0.95 

8 

6 

1.0 

2.024 

3.2824 

1395 

420 

2704 

403 

40 

33 

0.95 

16 

12 

1.0 

2.024 

3.2827 

5209 

1648 

10248 

2464 

41 

5 

0.98 

4 

3 

1.5 

1.996 

5.1580 

335 

92 

636 

97 

42 

34 

0.98 

4 

3 

1.0 

2.006 

5.1904 

449 

124 

856 

119 

43 

12 

0.98 

6 

4 

1.5 

2.006 

5.1939 

627 

182 

1206 

187 

44 

15 

0.98 

8 

5 

1.0 

2.010 

5.2105 

1388 

417 

2700 

459 

45 

35 

0.98 

8 

6 

1.0 

2.010 

5.2113 

1583 

478 

3082 

465 

46 

36 

0.98 

16 

12 

1.0 

2.010* 

5.2130 

6021 

1908 

11874 

2087 

47 

6 

0.99 

4 

3 

1.5 

1.950 

7.1511 

363 

100 

292 

104 

48 

58 

0.99 

4 

3 

1.0 

1.983 

7.2817 

491 

136 

940 

113 

49 

7 

0.99 

6 

4 

1.5 

1.984 

7.2923 

687 

200 

1326 

204 

50 

11 

0.99 

8 

5 

1.5 

1.997 

7.3459 

1113 

334 

2164 

356 

51 

59 

0.99 

8 

6 

1.0 

2.004 

7.3789 

1765 

534 

3446 

417 

52 

16 

0.99 

12 

8 

1.0 

2.004* 

7.3853 

3467 

1080 

6818 

1600 

53 

8 

0.999 

4 

3 

1.5 

1.319* 

17.2194 

475 

123 

916 

135 

54 

19 

0.999 

4 

3 

1.0 

1.484 

18.5617 

645 

180 

1248 

182 

55 

9 

0.999 

6 

4 

1.5 

1.365 

18.5919 

907 

266 

1766 

273 

56 

10 

0.999 

8 

5 

1.5 

1.635* 

19.9331 

1477 

446 

2892 

492 

57 

18 

0.999 

8 

5 

1.0 

1.807 

21.5183 

2047 

620 

4020 

618 

58 

21 

0.999 

8 

5 

0.5 

1.915* 

22.6040 

3401 

1032 

6688 

1200 

59 

22 

0.999 

8 

5 

0.3 

1.940* 

22.5819 

3843 

1162 

7520 

1400 

60 

17 

0.999 

12 

8 

1.0 

1.928* 

22.8096 

4569 

1428 

9022 

1876 

61 

20 

0.999 

12 

8 

0.75 

1.950* 

23.0445 

5235 

1636 

10332 

2295 

*For  these  NASTRAN  solutions  the  maximum  stress  did  not  occur  at  the  expected 
place  (o,R). 
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TABLE  3  -  PROBLEM  2  RESULTS 


a.  "SERENDIPITY"  ELEMENTS 


K 

KK 

N1 

N2 

a 

0  MAX 

SE  x  107 

GP 

NO 

DOF 

CP 

1 

5 

4 

3 

i 

17.345 

1.61776 

289 

78 

536 

76 

2 

9 

8 

6 

i 

23.853 

1.63167 

993 

296 

1904 

227 

3 

17 

16 

12 

i 

33.832 

1.64043 

3749 

1180 

7334 

1123 

b. 

TWO  OUARTER-PO I NT 

"SERENDIPITY 

'"  ELEMENTS* 

K 

KK 

N1 

N2 

a  °MAX 

SE  x  107 

GP 

NO 

DOF 

CP 

IS 

5 

4 

3 

1  73.936 

1.64582 

289 

78 

536 

76 

2S 

9 

8 

6 

1  94.836 

1.64769 

993 

296 

190 

227 

3S 

17 

16 

12 

1  131.71 

1.64869 

3711 

1168 

7260 

1106 

c. 

TWO  QUARTER- 

POINT  "SERENDIPITY 

"  ELEMENTS  AND  PERTURBED  TIP 

OF  CRACK 

K 

KK 

N1 

N2 

a 

aMAX 

SE  x  107 

GP 

«0 

DOF 

CP 

4S 

19 

16 

12 

1 

46.123 

1.71764 

3711 

1168 

7262 

1106 

*  A 

very 

minor 

change 

in 

the  algorithm 

for  a  in  the 

data 

generator 

added 

an 

extra  tier  of  12  elements  In  the  second  subdomain  for  the  third  solution  of 
Table  3b. 


to  elements  9  and  10.  Runs  3S  and  4S  gave  4.69723  as  the  value  of  from 
Equation  (8). 

The  stress  Intensity  factor  Kj  has  also  been  computed  by  fitting  the 
singular  behavior  of  the  solution  at  the  tip  of  the  crack.  Here  Equations  (1) 
and  (2)  were  used  to  obtain  Equations  (5)  and  (6)  which  can  be  used  for 
various  choices  of  grid  points  and  angles  6.  Figure  6  shows  the  numbering 
scheme  used  for  these  grid  points  in  the  computation  of  K^. 


Figure  6  -  Numbering  Scheme  for  Computation  of 
Table  4  gives  the  coordinates  (in  columns  headed  X  and  Y)  and  the 
displacements  (in  columns  headed  and  Uj)  of  the  grid  points  (in  column 
headed  PT)  required  to  compute  K^. 

The  results  from  Problem  2  are  presented  in  Tables  5  (all  "serendipity" 
elements)  and  6  (two  quarter-point  elements).  In  these  tables  the  column 
headings  I  and  0  refer  to  Fu(0)(I  -  1)  and  Fy(0)(I  -  2)  where  0  is  the  angle 
of  Figure  3.  The  grid  point  or  points  used  are  listed  under  the  PT  or  PTS 
heading,  depending  on  whether  the  one-point  formula,  Equation  (5),  or 
two-point  formula.  Equation  (7),  is  used.  For  the  two-point  formula  the  £ 
column  heading  gives  the  ratio  of  the  distance  of  the  closest  point  to  the 
furthest  point.  The  values  of  are  listed  under  the  heading. 


O'  Ul  O  W  ^  j  to  m  O'  U1  O  W  ^ 


TABLE  4  -  COORDINATES  AND  DISPLACEMENTS  REQUIRED 
TO  COMPUTE  THE  RESULTS  OF  PROBLEM  2 


a. 

PT  X 

ALL  "SERENDIPITY" 

Y 

ELEMENTS,  Nj 

U  x  108 

-  4,  N2  -  3 

V  x  107 

4  0.0000 

0.7750 

6.5910 

-2.2387 

3  0.0000 

0.8125 

4.0986 

-2.2318 

0  0.0000 

0.8500 

0.0000 

-2.2095 

5  0.0000 

0.8875 

0.0000 

-2.1496 

6  0.0000 

0.9250 

0.0000 

-2.0995 

1  1.0424 

0.8494 

1.8248 

-2.0111 

2  2.0847 

0.8488 

3.2685 

-1.8624 

b. 

r  X 

ALL  "SERENDIPITY" 

Y 

ELEMENTS,  N^  •  8, 

U  x  108 

Nz-6 

V  x  107 

0.0000 

0.8125 

4.6015 

-2.3021 

0.0000 

0.8313 

2.8135 

-2.2906 

0.0000 

0.8500 

0.0000 

-2.2710 

0.0000 

0.8688 

0.0000 

-2.2245 

0.0000 

0.8875 

0.0000 

-2.1784 

0.0265 

0.8498 

1.5205 

-2.1175 

0.0530 

0.8495 

2.6420 

-2.0149 

c. 

ALL  "SERENDIPITY" 

ELEMENTS,  Nj  -  16, 

N2  -  12 

PT 

X 

Y 

U  x  108 

r-. 

o 

X 

► 

4 

0.0000 

0.8313 

3.2419 

-2.3407 

3 

0.0000 

0.8406 

1.9829 

-2.3298 

0 

0.0000 

0.8500 

0.0000 

-2.3133 

5 

0.0000 

0.8594 

0.0000 

-2.2770 

6 

0.0000 

0.8688 

0.0000 

-2.2395 

1 

0.0138 

0.8499 

1.1071 

-2.2058 

2 

0.0276 

0.8499 

1.9128 

-2.1407 

TABLE  4  (Continued) 


d.  TWO  QUARTER-POINT" SERENDIPITY"  ELEMENTS,  N}  -  4,  N2  -  3 


PT 

X 

Y 

U  x  108 

V  x  107 

4 

0.0000 

0.7750 

7.2988 

-2.3260 

3 

0.0000 

0.8313 

3.4097 

-2.3460 

0 

0.0000 

0.8500 

0.0000 

-2.3551 

5 

0.0000 

0.8688 

0.0000 

-2.2649 

6 

0.0000 

0.9250 

0.0000 

-2.1944 

1 

0.0216 

0.8497 

1.7230 

-2.1774 

2 

0.0863 

0.8488 

3.7172 

-1.9273 

PT 

e.  TWO 

X 

OUARTER-POINT  " 

Y 

SERENDIPITY"  ELEMENTS, 

U  x  108 

«  8,  N2  =  6 

V  x  107 

4 

0.0000 

0.8125 

5.0238 

-2.3529 

3 

0.0000 

0.8406 

2.3685 

-2.3590 

0 

0.0000 

0.8500 

0.0000 

-2.3700 

5 

0.0000 

0.8594 

0.0000 

-2.2995 

6 

0.0000 

0.8875 

0.0000 

-2.2354 

1 

0.0135 

0.8499 

1.3565 

-2.2297 

2 

0.0540 

0.8495 

2.9315 

-2.0542 

PT 

f.  TWO-QUARTER-POINT  " SERENDIPITY 

X  Y 

’"  ELEMENTS,  Nj 

U  x  108 

=  16,  N2  *  12 

V  x  107 

4 

0.0000 

0.8313 

3.4982 

-2.3673 

3 

0.0000 

0.8453 

1.6639 

-2.3673 

0 

0.0000 

0.8500 

0.0000 

-2.3755 

5 

0.0000 

0.8547 

0.0000 

-2.3234 

6 

0.0000 

0.6688 

0.0000 

-2.2711 

1 

0.0070 

0.8500 

0.9697 

-2.2740 

2 

0.0281 

0.8499 

2.0931 

-2.1623 
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TABLE  5  -  PROBLEM  2  RESULTS  (ALL  "SERENDIPITY"  ELEMENTS) 


PT 

T 

a.  ONE  POINT  FORMULA  RESULTS  (5) 

. 

r  a 

1 

N1=4,N2=3 

V8,1^"6 

v16,v12 

3 

1 

3.14159 

3.86277 

3.97900 

3.85003 

4 

1 

3.14159 

4.46716 

4.52453 

4.45891 

1 

I 

1.57080 

3.67993 

3.49180 

3.71401 

2 

1 

1.57080 

4.52145 

4.42238 

4.53724 

1 

2 

1.57080 

3.71509 

3.79737 

3.60747 

2 

2 

1.57080 

4.38273 

4.69706 

4.09439 

5 

2 

0.00000 

1.82550 

1.66087 

2.01358 

6 

2 

0.00000 

2.56794 

2.15704 

2.89432 

b.  TWO  POINT  FORMULA  RESULTS  (7) 

PTS 

I 

e 

K 

V4»V3 

1 

Nl“8,N2*6 

^=16^2=12 

3,4 

1 

3.14159 

0.50000 

2.40364 

2.66199 

2.39939 

1,2 

1 

1.57080 

0.50000 

1.64831 

1.24517 

1.72656 

1,2 

2 

1.57080 

0.50000 

2.10324 

1.62532 

2.43196 

5,6 

2 

0.00000 

0.50000 

0.03308 

0.4630 

-.11251 
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TABLE  6  -  PROBLEM  2  RESULTS  (TWO  QUARTER-POINT  "SERENDIPITY"  ELEMENTS) 
a.  ONE  POINT  FORMULA  RESULTS  (5) 


PT 

I 

e 

N1=4,N2=3 

*1"8'M2,,6 

N1"16*H2-12 

3 

1 

3.13159 

4.68128 

4.59875 

4.56916 

4 

1 

3.14159 

5.81036 

4.87721 

4.80282 

1 

1 

1.57080 

4.62042 

4.59608 

4.55982 

2 

1 

1.57080 

4.98413 

4.96620 

4.92133 

1 

2 

1.57080 

4.76799 

4.81492 

4.77368 

2 

2 

1.57080 

5.73636 

5.34991 

5.01356 

«; 

2 

0.00000 

3.53653 

3.91476 

4.08655 

6 

2 

0.00000 

3.15213 

3.73316 

4.09673 

b.  TWO 

POINT  FORMULA 

RESULTS  (7) 

PTS 

I 

e 

S 

*1 

N1a,4,N2*3 

Nl“8,N23'6 

^-16,^2*12 

3.4 

1 

3.14159 

0.25000 

4.35219 

4.32029 

4.33553 

1,2 

1 

1.57080 

0.25000 

4.2567L 

4.22597 

4.19831 

1,2 

2 

1.57080 

0.25000 

3.79162 

4.27994 

4.53379 

5,6 

2 

0.00000 

0.25000 

3.92093 

4.09636 

4.07637 

5.  DESCRIPTION  OF  FEARS  COMPUTATION 


Because,  unlike  NASTRAN,  the  FEARS  program  has  a  built-in  data  generator, 
the  input  to  FEARS  consists  merely  of  the  basic  data  of  the  particular  problem 
to  be  solved: 

.  the  description  of  the  domain  and  the  type  of  boundary  conditions 
.  the  coefficients  of  the  pertinent  partial  differential  equation  of 
plane  elasticity 
.  the  type  of  error  norm  used 
.  the  type  and  amount  of  output  desired. 

The  error  norm  used  was  the  simplest  one  possible,  i.e.,  the  energy  error 
norm  L  for  q  -  2  and  q  *  °°.  The  input  parameter  p  set  to  1.0  specifies  the 

q 

L2  norm;  p  set  to  0.0  specifies  the  L*,  norm.  In  the  present  case  of  plane 
stress  the  strain  energy  density  is  defined  as 


where 


it 


a2  +  a2  -  2 ra  cr  +  2(l+y)a2 
xx  yy  xx  yy  xy_ 


a 

xx 


a 

yy 


r  3u  .  3vl 

L  ta  +  Y 

r^+  Y  3u] 
Uy  Y 


(9) 


(10a) 


(10b) 


o 

xy 


(10c) 


The  energy  norm  of  the  solution  s  * 
defined  by 


(u,v)  with  respect  to  the  norm  is 


E,2 


(11) 
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and  the  energy  norm  of  s  with  respect  to  the  L®  norm  by 


II-iIe,.-  II^IIl.  02 

T 

If  e  *  (ej,e2)  is  the  difference  between  the  exact  and  the  finite  element 
solution,  the  two  energy  error  norms  are  given  by  Equations  (11)  and  (12) 
respectively. 

FEARS  uses  elements  of  degree  1,  that  is  for  uniform  meshes  of  size  h  and 
a  smooth  solution  the  rate  of  convergence  in  the  energy  norm  is  0(h)  ■  O(^) 
where  N  is  the  number  of  degrees  of  freedom. 

6.  FEARS  RESULTS 

In  Problem  1  the  user  partitions  the  domain  ft  into  a  set  of  two- 
dimensional  subdomains  which  are  "curvilinear  rectangles".  The  particular 
choice  of  the  2-D  subdomains  can  influence  the  solution  in  some  way  because  it 
can  affect  the  aspect  ratios  of  the  elements  of  the  mesh.  This  partition  is 
characterized  by  the  coordinates  of  the  vertices,  the  curvature  of  the  lines 
joining  the  vertices,  and  the  numbering  of  these  vertices,  lines,  and  sub- 
domains.  Figure  7  shows  the  partition  of  the  domain  ft  used  for  r  *  .98.  The 
vertices  are  identified  by  numbers  in  circles,  the  lines  connecting  the 
vertices  by  ordinary  numbers ,  and  the  subdomains  by  numbers  in  squares .  The 
coordinates  of  the  vertices  are  given  in  Table  7. 

The  FEARS  program  then  constructs  adaptively  a  series  of  meshes  with 
respect  to  the  error  norm  selected.  The  FEARS  program  at  present  provides  the 
stresses  in  the  middle  of  the  elements,  and  the  stress  concentration  factor, 
SCF,  is  computed  according  to  the  choice  of  the  L2  (p-1.0)  and  L,.  (p=0.O) 
energy  error  norms,  respectively.  In  Table  8  the  data  for  the  next-to-last 
and  the  last  partitions  are  listed  under  the  column  headings  "FIRST  MESH"  and 
"SECOND  MESH”.  The  partition  described  under  column  heading  "THIRD  MESH"  was 
adaptively  constructed  but  not  computed  due  to  program  limitations. 

Table  9  presents  the  energy  error  norm  and  the  error  in  the  stress 
concentration  factor  (SCF)  for  solutions  computed  in  the  adaptive  mode,  the 
Lb  adaptive  mode,  and  the  solutions  computed  using  uniform  meshes.  Again  the 
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TABUS  7  -  COORDINATES  OF  THE  PARTITION  OF  DOMAIN  0  FOR  PROBLEM  1 


I 

X 

Y 

1 

0.0000 

1.0000 

2 

0.1000 

1.0000 

3 

0.2500 

1.0000 

4 

0.5000 

1.0000 

5 

1.0000 

1.0000 

6 

2.0000 

1.0000 

7 

3.5000 

1.0000 

R 

6.0000 

1.0000 

9 

6.0000 

0.0000 

10 

3.5000 

0.0000 

11 

0.9800 

0.0000 

12 

0.6930 

0.6930 

13 

0.4383 

0.8765 

14 

0.2377 

0.9507 

15 

0.0975 

0.9751 

16 

0.0000 

0.9800 

TABLE  8  -  HUMBER  OF  ELEMENTS  IN  THE  PARTITIONED  DOMAIN 
WITH  RESPECT  TO  BOTH  THE  ADAPTIVE  Lj  AND  L»  MODES 


Index  of 

2-D  Subdomain  Number  of  Elements  In  That  2-D  Subdomain 


FIRST  MESH 

SECOND  MESH 

THIRD  MESH* 

L2  L- 

L2  L» 

L 

2  oo 

i 

55 

148 

202 

154 

751 

559 

2 

52 

31 

139 

76 

511 

76 

3 

64 

37 

256 

73 

682 

103 

4 

67 

16 

238 

16 

841 

16 

5 

61 

16 

148 

16 

349 

16 

6 

31 

16 

67 

16 

202 

16 

7 

16 

16 

16 

16 

16 

16 

Total  number 

of  Elements 

346 

280 

1066 

367 

3352 

802 

DOF 

782 

636 

2242 

792 

*  Both  of  the 

meshes  under 

the  third 

mesh  heading 

were 

Adaptively 

constructed 

but  solutions 

mere  not 

computed  due 

to  program  limitations. 
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TABLE  9  -  ERRORS  IN  THE  ENERGY  NORM  AND  THE  STRESS  CONCENTRATION  FACTORS 


a.  FOR  MESHES  CONSTRUCTED  IN  THE  ADAPTIVE  L2  MODE 


NUMBER  OF 

DOF 

IWIe.2 

1  leME.2 

- 2 

ERROR 

SCF  ERROR 

ELEMENTS 

x  104 

1WUe,2 

x  100 

IN  SCF 

SCF 

x  100 

346 

782 

1.1119 

15.39 

0.0922 

4.58Z 

1066 

2242 

0.6476 

8 .96% 

0.0325 

1 .62? 

b.  FOR 

MESHES  CONSTRUCTED 

IN  THE  ADAPTIVE  L»  MODE 

NUMBER  OF 

DOF 

11*11 E, 2 

IWI..! 

ERROR 

SCF  ERROR 

ELEMENTS 

I1UIIE  2 

IN  SCF 

SCF 

x  lo4 

x  100 

x  100 

136 

342 

1.8876 

26.14 

0.1625 

8.08Z 

280 

636 

1.6874 

23.36 

0.0931 

4.63? 

367 

792 

1.5875 

21.98 

0.0501 

2.49X 

c.  FOR 

UNIFORM  MESHES 

NUMBER  OF 
ELEMENTS 

DOF 

IMIm  . 

IWI.,2 

T^C 

x  100 

ERROR 

IN  SCF 

SCF  ERROR 

SCF 

x  104 

x  100 

112 

290 

1.9647 

27.21 Y 

0.3208 

15.96Z 

448 

1026 

1.0866 

15.05? 

0.1077 

5.36Z 
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•treM  concentration  factors  are  computed  by  extrapolation  from  atreaaes  in 
the  center  of  the  elements. 

The  energy  error  norm  ||  e  ||  ^  can  be  easily  computed  if  the  exact  energy 
is  known.  If  is  the  exact  energy  and  e  is  the  energy  of  the  Finite 

Element  solution,  then 

He^E,2  "  ®EXACT  "  6 

When  the  NASTRAN  and  FEARS  results  are  extrapolated  to  the  limit  ~ 

-7  EXACT 

«J>. 2130150  x  10  .  Similarly  by  extrapolation  to  the  limit  we  obtain 

SCI?EXACT  "  2*0106.  (The  tilde  notation  indicates  the  extrapolated  values  of 
eEXACT  and  SCFEXACT‘) 

Table  10  gives  the  error  estimators  for  the  energy  norm  |e|  with 

E 

respect  to  the  adaptive  mode  and  for  uniform  meshes. 

Table  11  gives  the  error  estimators,  scaled  in  the  same  way  as  the  stress 
concentration  factors,  for  the  stress  concentration  factors  with  respect  to 
the  L„  norm  and  for  the  stress  concentration  factors  for  uniform  meshes. 

For  consistent  comparison  the  error  extlmator  for  the  energy  norm  must 
be  compared  with  the  actual  error  measured  in  the  norm.  Table  12  shows  the 
actual  errors  in  the  energy  at  the  point  (0,.98)  and  the  error  estimators  with 
respect  to  the  adaptive  mesh  and  for  uniform  meshes,  although  the  maximal 
error  has  not  necessarily  occurred  at  this  point. 

Figure  8  shows  a  partition  of  the  domain  fl  for  Problem  2.  As  in  Problem 
1,  the  vertices  are  identified  by  numbers  in  circles,  the  lines  by  ordinary 
numbers, and  the  subdomains  by  numbers  in  squares.  The  coordinates  of  the 
vertices  are  given  in  Table  13.  The  FEARS  program  requires  that  a  vertex  be 
placed  at  the  top  of  the  crack. 

Because  the  stresses  are  infinite  at  the  tip  of  the  crack,  the  only 
adaptive  mode  that  could  be  used  is  that  of  the  L2  energy  norm.  By 
extrapolation  as  before  the  exact  energy  of  the  solution  was  found  to  be 
1.649701  x  10  .  The  energy  release  procedure  for  computing  the  stress 

Intensity  factor  Kj,  using  Equation  (8)  with  respect  to  the  norm,  requires 
adaptivity  with  respect  to  that  norm. 


TABLE  10  -  ERROR  ESTIMATORS  AND  THE  ACTUAL  ERROR  IN  THE  ENERGY  NORM 
a.  FOR  MESHES  CONSTRUCTED  IN  THE  ADAPTIVE  L2  MODE 


fkimber  of  Elements 


DOF 


'  "E.2 
x  10*  • 


Estimator 
x  10* 


Estimator 

TUTT^ 


346 

782 

1.1119 

1.0657 

0.9584 

1066 

2242 

0.6476 

0.6461 

0.9977 

b.  FOR  UNIFORM  MESHES 

Number 

of  Elements 

DOF 

"•"w 

x  10* 

Estimator 
x  10* 

1  a  _  Estimator 

112 

290 

1.9647 

1.6878 

0.8591 

448 

1026 

1.0866 

1.0293 

0.9472 

TABLE  11  - 

ERROR  ESTIMATORS  FOR  THE 

STRESS  CONCENTRATION  FACTOR 

AND  THE  ACTUAL  ERROR 

a . 

FOR  MESHES 

CONSTRUCTED  IN 

THE  ADAPTIVE  L„ 

MODE 

Number 

of  Elements 

DOF 

Error  in  SCF 

Estimator 

Estimator 

Error  In  SCF 

136 

342 

0.1623 

0.3251 

2.0006 

280 

636 

0.0931 

0.2090 

2.2448 

367 

792 

0.0502 

0.1886 

2.6617 

b. 

FOR  UNIFORM  MESHES 

Number 

of  Elements 

DOF 

Error  In  SCF 

Estimator 

Estimator 

Error  in  SCF 

112 

290 

0.3208 

0.5251 

1.6369 

448 

1026 

0.1077 

0.3462 

3.2144 
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TABLE  12  -  MAXIMAL  ERRORS  IN  THE  ENERGY  AND  THE  ERROR  ESTIMATES 


a.  FOR  MESHES  CONSTRUCTED  IN  THE  ADAPTIVE  L ®  MODE 


Number  of 
Elements 

DOF 

3 

Actual  Error  x  10 

3 

Estimator  x  10 

Estimator 

II*  lE.» 

136 

342 

2.2618 

2.0987 

0.9279 

280 

636 

1.2365 

1.3493 

1.0913 

367 

792 

1.1976 

1.2172 

1.0164 

b. 

FOR  UNIFORM  MESHES  WITH 

RESPECT  TO  Loo 

Number  of 

DOF 

3 

Actual  Error  x  10 

3 

Estimator  x  10 

o  Estimator 
6  IUM 

Elements 

1 lel Ie,« 

112 

290 

3.9371 

3.3896 

0.8609 

A  48 

1026 

2.4662 

2.2349 

0.9062 

TABLE  13  -  COORDINATES 

OF  THE  PARTITION 

OF  THE  DOMAIN  0  FOR  PROBLEM  2 

I 

X 

Y 

1 

0.0000 

1.0000 

2 

0.2000 

1.0000 

3 

0.8000 

1.0000 

4 

2.0000 

1.0000 

5 

6.0000 

1.0000 

6 

6.0000 

0.0000 

7 

0.7000 

0.0000 

8 

0.4373 

0.5466 

9 

0.3131 

0.6261 

10 

0.1373 

0.6864 

11 

0.0000 

0.7000 

12 

0.0000 

0.8500 

13 

0.1667 

0.8335 

14 

0.3801 

0.7603 

i 


I 
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TABLE  14  -  ERROR  ESTIMATORS  FOR  MESHES  CONSTRUCTEE  IN  THE  ADAPTIVE  L_  MODE 


Step 

Number  of  Elements  in  2-D 

's 

Number  of 
Elements 

1 

2 

3 

4 

5 

6 

7 

1 

16 

16 

16 

16 

16 

16 

16 

112 

2 

16 

19 

16 

16 

16 

19 

19 

121 

3 

19 

19 

16 

16 

16 

19 

19 

124 

4 

22 

25 

16 

16 

16 

43 

34 

172 

5 

25 

25 

16 

16 

16 

45 

34 

175 

6 

40 

37 

34 

16 

31 

64 

46 

268 

7 

43 

37 

34 

16 

31 

64 

46 

271 

8 

55 

46 

46 

16 

43 

67 

58 

331 

9 

58 

46 

46 

16 

43 

67 

58 

334 

10 

187 

112 

160 

16 

145 

238 

166 

1024 

11 

187 

118 

160 

16 

145 

238 

166 

1030 

12 

190 

118 

160 

16 

145 

238 

166 

1033 

DOF 

£  X  10^ 

HellE,2 

4 

x  10 

'  6  S,2  s 

Estimator 

lluflE,2 

x  100 

274 

0.1569 

0.8975 

22.09Z 

0.6881 

286 

0.1579 

0.8405 

20.69Z 

0.8019 

292 

0.1583 

0.8152 

20.07Z 

0.7953 

376 

0.1604 

0.6796 

16.73Z 

0.8776 

382 

0.1606 

0.6631 

16.33Z 

0.8756 

566 

0.1623 

0.5177 

12.74Z 

0.9323 

572 

0.1624 

0.5063 

12.46Z 

0.9350 

674 

0.1628 

0.4611 

11.35Z 

0.9473 

680 

0.1629 

0.4547 

11.19Z 

0.9493 

2024 

0.1643 

0.2647 

6.51Z 

0.9626 

2034 

0.1643 

0.2599 

6.40Z 

0.9878 

2040 

0.1643 

0.2568 

6.33Z 

0.9900 

Several  versions  of  Problem  2  were  solved.  In  the  basic  case,  the  tip  of 
the  crack  was  placed  on  the  nodal  point  with  coordinates  (0.,.85).  Other 
cases  were  obtained  by  varying  the  y-coordinate  of  the  tip  of  the  crack  from 
.85.  (In  the  FEARS  input  this  perturbation  is  effected  by  changing  only  one 
number . ) 

The  presence  of  the  singularity  causes  a  special  situation  in  the 
following  sense.  The  construction  of  a  completely  optimal  mesh  sometimes 
involves  the  refinement  of  a  very  small  number  of  elements  (possible  only  one) 
at  one  step  of  the  FEARS  procedure.  This  situation  could  be  avoided  by  using 
another  FEARS  command  which  increases  the  number  of  elements  to  be  refined 
essentially  without  additional  computer  cost. 

Table  14  shows  how  the  sequence  of  meshes  was  created  when  complete 
optimality  of  the  meshes  was  desired.  Table  15  gives  the  basic  results  for 
the  Finite  Element  solutions  of  Table  14. 

Tables  15a  and  15b  show  other  sequences  of  meshes  adaptively  constructed 
by  strategies  which  avoid  a  small  increase  in  the  number  of  degrees  of  freedom 
during  mesh  refinement.  These  strategies  consist  of  inserting  two  (Table  15a) 
or  three  (Table  15b)  so-called  "short  passes"  as  required  into  the  usual 
adaptive  mesh  refinement  procedure.  Table  16  shows  the  performance  of  the 
uniform  meshes - 

The  stress  intensity  factor  is  computed  by  the  energy  method  which 
requires  the  meshes  to  be  the  same,  topologically  speaking.  The  special 
features  of  the  FEARS  program  provide  a  means  for  obtaining  consistent  meshes. 
The  length  of  the  crack  can  be  increased  very  slightly  by  altering  the 
y-coordinate  of  the  vertex  at  the  tip  of  the  crack.  The  meshes  can  be  so 
constructed  in  an  adaptive  manner  as  to  produce  the  same  mesh  (topologi¬ 
cally  speaking)  as  before.  The  number  of  the  adaptive  steps  used  to  con¬ 
struct  this  second  mesh  depends  on  how  much  the  crack  length  is  increased.  A 
very  small  increase  in  crack  length  extends  the  number  of  admissible  steps  and 
decreases  the  error  caused  by  approximating  the  derivative  by  a  difference, 
but  it  also  could  increase  the  effect  of  round-off  error. 

Table  17  shows  the  maximal  number  of  elements  in  the  topologically 
identical  meshes  for  the  basic  case,  the  energy  of  the  Finite  Eler  t 
solution,  and  the  growth  of  the  error  estimators. 
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TABLE  15  -  MESHES  CONSTRUCTED  IN  THE  ADAPTIVE  h  MODE  WITH 
CURTAILMENT  OF  SMALL  REFINEMENTS  1 


a.  BY  USE  OF  TWO  SHORT  PASSES  BETWEEN  LONG  PASSES 
Step  Number  of  Elements  in  2-D's  Number  of  DOF  e  x  10^  1 ! e I i e  2 


I 1*1 lg  ?  _  , 

tt-  I  r»  -  9  -  Eat  1— tour 
'  E,2  Tlej  / 


1 

2 

3 

4 

5 

6 

7 

Elements 

x  104 

1  ’  E,2 

x  100 

1  1*11 

1 

46 

46 

46 

16 

46 

64 

46 

310 

640 

0.1616 

0.5841 

14.38% 

0.7806 

2 

49 

58 

46 

16 

46 

64 

55 

334 

688 

0.1624 

0.5079 

12.50% 

0.8979 

3 

70 

73 

55 

16 

52 

91 

79 

445 

880 

0.1632 

0.4190 

10.31% 

0.9357 

b.  BY  USE  OF  THREE  SHORT  PASSES  BETWEEN  LONG  PASSES 


Step  Number  of  Elemnts  in  2-D’s  Number  of  °°F  e  x  10®  llell 

Elements  L 

1234567  *  10 


11  *E,2 
x  100 


Estimator 

TlelL'  ' 


1  46  34  16  16  25  55  37 

2  157  151  97  22  112  223  157 


474  0.1613 

1740  0.1640 


0.6094 

0.3099 


15.00%  0.9019 

7.63%  0.9713 


TABLE  16  -  ERRORS  AND  ERROR  ESTIMATORS  FOR  THE  UNIFORM  MESHES 


Step  Number  of  Elements  in  2-D’s  Nv“ber  of  DOF  e  x  10 
-  Elements 


I -I  I  11  1  ^  Estimator 

"y  TMT^e‘TRI^- 

104  x  100 


1  16  16  16  16  16  16  16 

2  64  64  64  64  64  64  64 


0.1569  0.8975  22.09*  0.6881 

0.1616  0.5809  14.29%  0.7197 


TABLE  17  -  ENERGY  AND  ESTIMATED  ERROR  IN  THE  ENERGY  WITH  REFERENCE  TO  THE 
LOCATION  OF  THE  TIP  OF  THE  CRACK  FOR  TOPOLOGICALLY  IDENTICAL  MESHES. 


Y 

Number  of 
Elements 

DOF 

£  x  106 

Error  Estimator  of  the 

Error  in  the  Energy  x  10° 

0.8600 

112 

274 

0.1601 

0.4131 

0.8400 

124 

292 

0.1552 

0.3900 

0.8450 

124 

292 

0.1567 

0.4049 

0.4860 

124 

292 

0. 1571 

0.4079 

0.8540 

172 

376 

0.1617 

0.3675 

0.8492 

271 

572 

0.1621 

0.2224 

0.8501 

1030 

2034 

0.1643 

0.0660 

TABLE  18  - 

Number  of  Elements 

ENERGY  AND  COMPUTED  ESTIMATE 
FOR  THE  CASE  OF  y  =  .85 

DOF  e  x  106 

OF  THE  ERROR  IN  THE  ENERGY 

Error  Estimator  of  the  g 

Error  in  the  Energy  x  10 

112 

274 

0.1569 

0.3814 

124 

292 

0.1583 

0.4206 

172 

376 

0.1604 

0.3557 

271 

572 

0.1624 

0.2240 

1030 

2034 

0.1643 

0.0659 

TABLE  19  -  STRESS  CONCENTRATION  FACTOR  K  COMPUTED  BY  THE  ENERGY 
RELEASE  APPROACH 


Y 

Number  of 
Elements 

DOF 

K1 

Estimator  for 
the  Error  in  Kj 

Estimator 

Y  x 
K1 

0.8600 

112 

274 

4.3928 

0.2116 

4.81% 

0.8400 

124 

292 

4.3036 

0.2082 

4.83% 

0.8450 

124 

292 

4.3479 

0.2116 

4 . 86% 

0.8460 

127 

292 

4  3569 

0.2122 

4.87  % 

0.8540 

172 

376 

4.5175 

0.1916 

4.24% 

0.8492 

271 

572 

4.5961 

0.1337 

2.90% 

0.8501 

1030 

2034 

4.7224 

0.0365 

.77% 

34 


Table  18  gives  the  energy  and  the  error  in  the  energy  for  the  meshes  of 
Table  17  for  the  basic  case  of  y  »  .85. 

Equation  (8)  is  now  used  to  compute  the  stress  intensity  factor  Kj  by 
replacing  the  derivative  with  the  one-sided  difference.  Table  19  gives  the  Kj 
factors  and  the  computed  estimators  for  various  values  of  y. 
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-  O  C: 


APPENDIX  -  PROGRAM  LISTING  OF  DATA  GENERATOR 


PRCG-AM  J AT  Ail  l INPUT , OoT°UT , T  A°E5* INPUT  » TAPEfexCUTFUT ,TAPcl«TAPE2, 
1  TMaE3,TAPE'«,TAPc5) 

ClMENSICN  rrai 

HOLE  PR02CE1  DATA  GENERATOR  FOR  NASTRAN 

T «PEl=iR ID  SARDS 
T A°c2sElfc  mENT  (CIS’So)  CARDS 

TAP£3=;.NS T-AJNT  cas:s 

TAO£i*=FJi?Ci  CLPSS 

CIME N$  I CN»N2* 1 

iOh  MoA*  3£U(31)  ,  SI0£Y»3i) 

DIMENSION  >N 3  (SEE  FORMAT  NO.  fc°fc) 

N3  =  ND M-Jca  OF  ELEMENTS  IN  X-DIRECTICN  IN  TMIKC  sUS-COHIN 
CCi-MCr.  *1 OT  J  .,) 

CIM£NSI0a»2*N2»1  (SEE  FORMAT  NO.  25m) 

'.O.-mCn  iAAlol  )  « ij  J  t  b 1 )  ,P(61) 

OMENSIONM  (SEE  FO-MAT  NC.25TT 
i.s  (  HUM  l£R  CF  cLi.Mci.rS  A  »  C  L  N  _  akC  iri  FIRST  SU3“v/CM  LIN)  ♦  1 
wOKMON  TOP*  <i.„.),TM£T  A  ,XL3( AS-.  >,SCU-S> 

G-n  POINT  -RRAYS 

common  xiioaoo)  .YUoaoo) .NENiiaoao) 

ARRAY  DIMENSIONS  TO  3E  CHECKED 
JmT  a  NuE.NN.MF.NLIMJ/ 30, 1030,10000, 2000/ 

INFUT  DATA 

MsNUM3Ek  OF  ELEMENTS  ALONG  Y-AXIS  OF  FIRST  SU3-D0.1AIN 
N2  =  MJMoER  OF  ELENtMS  IN  Y-DIRECTION  IN  THIRD  SUb-D CHAIN 
MsPASIUS  OF  CIRCuE 
HsmcI^MT  OF  RIGHT  SIDE 
Xl=wcNGTH  OF  1CP  SICE 

AMsmSPECT  ratio 

-'EL0(S,*IN1,N2,A,H,XL,AR,KK 
IF(K«i.LT.llKK*l 
nnl  T'.  (o  •  *9)  HI  •  N  2,R  «H«XL.  AR  «KK 
Ac  FCM<-T  i  1X,*N*  ,N2,-\»M,XL,AP,KKS*,tU3t^E15.E,I10) 


ChESK  N2 

IF  (N2.LE.NN2)  GO  TC  255 
nn.1  TS  (a  ,  25AJ 

25-  FC; MAT (1X,*N2  TCO  DIG*) 
GO  TO  993 

253  IPI0*1 
IEI8«0 
IDl»3 

isioi*u 

cosi*i.o 

ISI02«21 

wIH2«t*N2»l 


*2' 


<^v 

IT* 

& 


FHlCki)h>a  PAuS 


X2t.l***Ml 

*m=ru 

XN2=N2 
Nl»C  =  'U*l 
N2J-C i=N2»l 
LlK12=2*N14l 
X2N2=2*N2 

^1=3,141532 65 35a 97 93 

n:;.~U  *0-  „LFMA 

anGIt=<°I/U0.0)*27.S 
.'C  A.'J  U  1=1  *  1.  V  te 
5TGS£=A**<xfU/XN2)*ANGIT4l.O 
3TGeL  =  SJC9£«<i7/Hl 
ANGNE H= AT AN2< 1.0, STOKE! 

J  =  1 

IE  (H^  (  Uf.GI  T-ANo  4 £.U  1  /A'iGITI  .vt..  0011  uO  TO  152 

150  «NoIT=ANGNEW 
«F  i  T  E  ( o  « 151 ) 

151  FO-MAT t lHi/4; X,*ALFHA  013  NOT  CCNVi *G: t  !  Itt  !  <*! 

4U  T  J  £55 

l'G^  «c.Tr:  (i,15otir.C-NGM,  J 

15  1  H)s.1i.T(J]X,*iLfHU,i(C  IT  E« ATIQNS  -»  ,  •  lo  .  «  ,  1101 
DETERMINE  SUE  CC?*4XN  X 
GAMK A=ANGNFM 

CIO  T=  IXN1/XN2)  *K*GAK)<A*AK 
?'j:ST  =  ?.OIST 
tLFHA=ATAN2IH,R*0ISTI 
XL  Zc  (<0-H”R 

XLt  =?CiST/C  JS  ( ALPHA  I  -F. 

XL.  If  F*XL.,-XlZtfiO 

lm.TAH=.5»PI-ALPHA 

fHi.TA«ll=a.O 

XL  S ( 1 !  =  XLZE3C 

LO  -i?  1=2, 100 

wl*l*  I  XLSU-U  /XNl)  *Aft 

1HL  Trill  I  «SIII  /  F  ♦  T  He  T  A 1 1  *1 1 

N=1 

XF  WMiTAII)  .LT.THiTAMl  GO  TO  436 
GO  TO  A o i 

436  XX  =  TH;1Am/C  THETAh*TANI  7A.FHA)  1 

XlSCI  =  SGr.7<<xX-.-*iIMTHETAm  )>  **2  «  lh-4*CCS  ITHET  MI  II  1**2! 
A3?  CONTINUE 


-.3?  ;F  IN.lE.NM  GO  TO  256 
MAiTt  16,257) 

257  FOf.HAT  11X,*T0PX,  THETA,  XLS,  ANC  S  CIPENSION  TOO  SHALL*! 
GO  TO  5^9 


256  GCALtF=rHETAM/ThETACK) 

UO  439  1  =  1, N 

ri39  THETA<:)sSCALEF»THETA«I) 

tofxn=«oist 

00  <*40  1*1. H 


o*2**x 


38 


ooo 


<*4  0  TCPX(l>*tT«eT  AU>/T*-ETAm*TOPXN 


SU6C<3HAI<I<  I 


IP=0 

NEWIP=0 

00  5  -  i.  L  =  1  •  N 

XIF=  jlNlTHET*.<L»»*R 

XX^=:C3<Tn<TA<L>)*3 

AlNC*<T3oxiD  -Xl?)/XiMl 

YIi«C*»H-YXP»  /X2K1 

GC  *59  I=1,-IH12 

XP=IF*1 

NEWI°*NE AIP+1 

NEMiPI snEMIP 

X<iP)=X  JF+FLOA  MM  I  *XINC 

Y<.r-)  =  <IP*f-C6T<i*n»YXNv 

IF  <<<•  GT  .1  .ANC.  ic.t-O.K<-l)  GC  30  323 

GO  TO  324 

323  HiP>s''tlPt*Vl*C/Z* 

GO  TO  49s 

324  IF  (KK.GT.l  .ANO.  IP.tC.KKUl  GO  TO  325 
GC  TO  493 

325  IF  GO  TC  499 

»i:B|*YtiP»->IAC/2. 

459  -All  W5T“Tl<i,XF.X<IF>*Y(IP>> 

X£Kia=NEHlP*MPO 

SOU  CONTINUE 
NEHP*LiM12 
LC  G50  L=2,N 

-NGLi=i.5*<TMtTA(L> ♦TH£TAIL-1»» 

IF=IP«1 
NEKIF=A£  AIF41 
>*£*>  tiPI  *NtNIP 
XUP»*SISiANGL6>*R 
vap»*cos(ANGt.ei*ft 
1ST  0H£* IF 

-ALL  N5TFTU1.1P.X(IP>,Y«IPII 

00  >25  1=3*LIH12*2 

IP=IP*1 

KE4lP»NtYlP*l 

Nto  <  XF • *NEHIP 

IA=  t‘.-2>  *LMl2*l 

19= I L -l J  *LlNl2 ♦  I 

XUP>*0.5MX(IA»*XIIBH 

Y«IFI *Q,3*(Y(IAJ*Y(IQI » 

IF  IKK.GX.l  .ANC.  w.ELi.2  .AND.  X.tQ.NlFOl  GO  TO  360 
uC  to  351 

350  X(ie«s(3.0*XllA)«X(lB>l/4. 
r<;P)a(3.-»Y<IA)»Y(ie»l/4. 

351  -ALL  HRTPTUI  *1°*XXIP»  »fUPI» 

XGHIA 

IE20«ieiJ«l 

IG2*XA*2 

AG3*I6'2 

164*16 

aG>*IA*1 
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)  t » 


IGb*lSTO>(E 

lG7=ia-i 

1G0-IP 

•i-.iT-  <2,^321*6  lC,Ir-4l,»jl,Ik>2,XtaJ«IC->t*G5t  IGto, I EXC, I £IG»  XG7  t  IGE 1 
1  JJ1 

PCf  1M  (4H.:S2C»  'oI?.l«», 17/144, 17,3X0 

->2E  ISIO^tsIP 

kEkI-=N£HF4wIM2 
•j  5  U  *-  c  i  •  T  1 N  0  £ 

^  T  a  X  *  I  P 

p‘zf*zl  fuj  cf  sp(.i  ;ari3$ 

Xl  =1 

h-.XTE  ( 3,30  0  «  USIOl.IC.IFP) , 1PF=*K ,LIN 121 
303  FO-MM  UHSFCl«<tX,4I8l 

XOaJGPaXN  II 


PNGc  £»Al-HP/XN2 
»If.C  =  H/XN2 

::  39r.  x=i,  m2  pc 

:•£!  a  ( 1)  =FL3AH  I-l)#fiFGL£«rHETA« 
i'z'.  ill-'T  <  i  >  =H-FuCPT  (1-1) *YI  iC 
NCfcIF=NP*ai*MrC 
IC  iJj  L=2,N2FC 
XIP  =  aINC:EHU»»*^ 

YIr  =  CCS  l=£TA(LI>*S 

*IKC*(R0:5T-»IPl/x2M 

YI62=  (aXOEYUJ  -YIPI/X241 

»F  (£.aa.N2Pu)  GO  IC  597 

YIP=0.Q 

YINC*0.0 

597  GO  595  I  =  1,».IH12 
XF  =  X  F* 1 
NEwiP=f4EaiPn 
rt£ • ( X  P I =  N£  KIP 
X (IP  I  =  X IP 4Fi.C a  I  a- 1»*X INC 
YUP*  =YlF»PLOAia-l>*YIN0 

sS9  :u.l  wptptkx  ,xp,x(xp»  ,vcif»» 
N£^I  =  =  N£XXP*.41P0 
60  u  ONT  IN  J  c 

iAGI'l  =  Cn-il*LIf1i2 
:.£wip=,\prsx 
20  625  u=2,N2PC 
ANit£*0.5*C3t7  4  «L»*BETAXw-l»» 
IP=IP«1 

xCXP»=5XMPNGtH>*fi 

KEHF  =  KGHC»1 

i.£t»liPM«ttfIP 

Y(XFI=C03(ANGLE)*R 

call  K5T’iia,iP,xaFi»Y«XPH 

XSl OPE  =  I? 

I8i.T.<T  =  KPTS>I»  a-2»*LlrtlZ 
00  62«,  I  =  3,LXM12,2 
IP= IF ♦ 1 
KE*IP  =  NE*4IP*1 
NFMXFI  *N£WIP 
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U«i»ST«T*I 

*.8=I^ST*T»3 

CALL  MrT'»THI»1P»***f,»,<XPI> 

1.21*1* 

lG2=XA-2 

:g*-*I2 

;gv*ia-i 

*G£=IST03£ 

1  101 

el*  ISIO-ve.1!? 

fc25  iAi.TrT*I  iSf<<T 

SG3C3«AiN  HI 

i.fTSlI*!? 

;.=xu-S3ISl 

U3*0/«<H/XN2>*AR» 

XN3*N3 

Hp3*A3*1 

IP  <LlP3.U.4tIM3>  GC  TO  G« 

o96  «««*"  100  S"LL'’ 
ttO  TO  999 

e95  XSgT (II *AOXST 
XlNO*C/XN3 
NE«tf:*NPT31I*N2«> 

S^AOISt^toJuL-lJ’XlNC 
X3CT  Ct»**IF 

Ylk  *0  •  0 
YlWC*H/X2Ni 
CO  699  I=1.LIH2 
I  P*Ip*l 
KE>iIP*t*E*tXp*1 
NiKUPi*NE«P 

?55»*5xwc*i«-n*rxHC 

-cc  1. .  B^TPTHs9S»Xc»X**r‘,<,IP" 

KeiiP*Ne«iMN2fo 
700  CONTINUE.  T11,, 

1ST  OVE*N«>TtI*  NJ*UIH12 
00  *710  I*1*GIP2»Z 
lftfc<II*ISTOR£ 

710  iSTOf*i*XSTCRt.-CUl< 

HAli-luZl***1-1*12 
ISJw*?E*NPTSlI 
00  Eli  I«2»UNi»2 
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1AA ( X ) =  1  STORE 
711  i'TC;  ^  =  ljTjRi-MPO 
£  t  **r=NrI5*I 

7iJ  l=2,_IFl 

IP=IP ♦! 

FEWlP  =  U£.(IP*l 

:.£  M  Ir)  =  N£»Ir 

iKV.si  .5*  (XECT  (i.1  ♦X3UT  <(_-!)  ) 

x(i >i =iro?£ 

y  <!  =  j  =  j .  j 

v A OMPTl  (:>“,I=,  X(1F>  ,Y  (IMI 
LC  77,  islitlHX 

77C  13ji:>-NPT$Xi*<L-2l*LIM2»I 
20  775  I* J» tiM2 »2 
XP=I=*1 
NXAIPsilmXF*’! 

'.;K|  1 r )  =  N i  W I F 
14=144(11 
iJ=Ic=( I) 

*l*3l=GTC?t 
r  (!*-'»  =y  ( ld» 

'i.L  w;T3Ti(  j9=,1P,X(IP>  ,Y  (I?l  ) 

*01-14 

*  £  *  0=  xr.  I  £♦  1 

iG*  -  1 4  4 ( X-2) 
lG4  =  i?ri 1 1-2) 

IG**  =  I6 

iG5»IAAU-ll 

•  Go  -  *  S T  Q  ■<  £ 
iG/^iaa (i-ii 
XGt=IP 

"FiTt  (2.2u2IIEIC  ,IPIC,TG1,  IG2,  XCX, IGA, 135,  1G6, I EIC.I E J C, IG7 , IG 8, 
1  101 
775  IS)  0*E  =  I F 

-0  7-*9  J  =  1«LIK2 
7 1*5  1 4 A  (  J 1  =  I  J 3<  J) 

t'E«IP='(£i(IP'*LIt,2 
750  CCI  TIMJE 

ChCCK  OkXu  PCXM  AFRAYS 
*F  (Ir.CC.NI PI  GO  10  7**3 
W-ei  TE  <5*  r*7) 

7£i  FC?-I- T(IX**X,  Y  ClKEFilONS  TCC  SMAU.»> 

GO  TO 

p'tP4->’£  REST  CF  5PC1  C  AaDS 
7*.-  IC  =  2 

40  735  IQ=HP12  *ISTOS€ 

IF  (((1*1  .ME. 0.0)  GO  TO  785 
*.PH£(j.jJo)IS1C1«IC*I0 

785  COM  I NUE 

F' £  F  A  R  t  FOKCE  CAROS 
CO  736  X -1 • wXP2 *2 

786  F«I)=2.o 


fvnC  d°** 
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Fli»*l.l 

FUl»*2>*l««  .  . 

OC  ?■>?  I*2»tlfl2«2 

li  7 

rh-  t*-b  •  i*XN2/H 
OC  7o%  i3ltLlM2 
7oc  F  (il  =  F  1 /FNkM 

ISTOJcsNp^-‘**^*L**^ 

.i.  F  =  iSTuRi-cI<'2*l 

CO  Tea  1U® » ISTOSe 

roe  -iriTtlH.wiOISlCZfiOtF*!’ [fj?/*F8.3> 

^U«.  FOf  MATOHFOR-i  *21 

j  n0^4(.  T£PHNMIOk  CF  01T*  6ENE«*TC» 

IF  (ir.t-C.0--U>  GO  1C  l'5-‘ 

„,.,Y  CftU  POINTS  t»*aao>*> 

3HJ  FCiKMATUU/H3X»*TOO 

oC  TO 

150  CONTINUE 

ta.olillF  .It-IC 

Wr'iTtU.l7S><l.N  =  *‘l‘jI*i*IP‘ 

,  o-  t.»».  c»»i 

DO  105  l*t«*fcK 
FfcFOU*  lui»»H  ,t«. 

FOfHAl  C2AX»6H  ,8X/3X*2IS 
LO  102  1=1*5 


10c 

105 

110 


J=MK» 

HIKt=NEMiJI 

WniTi<a.U0»H 
FC-NAT (5HC132C8 
PtNINO  1 

r-  EfclNO  2 
AENIND  3 
•.ENlNO  <* 

9enino  a 

STOP 


*lbX»6ia/1H**TX,2*8> 


c 

C 

C 


ENSURE  RUN  NllP  AoCRT  RtClUSE  C*T«  IS  -»C 

<j5J  RENlNC  1 
BEHIND  2 

F.EvUNO  3 
REnINO  •* 

ENi-FI-E  1 

CNuFIEE  2 

iU  *<■>-««•• 
STOP 


END 
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iU‘90uTI:t£  W?TFTl<i,IP»X*Y» 


K^lhS  A  G«IC  ;ARU  OK  TAFtl  USING  OKE  CF  THREE  FC’MATS 

irca=j 

*F  (1  ijTtl)  3  0  TO  **41 
«=1 

c  =  ;  T<  X»X*f»Y> 
fi-J.  1  •«  1 1  i  ^  o  S  3  -)7  2  2 
.  =  (  1*  J  ./  3I)  *6  T  AN2  <  Y  ,X  J 
-RITE  U,  Unix?.  l„Or,F,C.  IP. 1? 
i  u  1  ECr  *4  I  (  it  riv>~  I  D  »  **  4  »  I  ^  f  IX  «  X  f  <«X  •  F  $  •  7  t  F4  •  4) 

441  .F  (X.jE.1.1  GC  TO  442 

»41TEU»li4)IF,XC0l.,X,Y,IP,tP 
IOl  r  0-  MAT  i  4HGHI j  ,<.X,Z?,«X,F?.7,F8.6) 

i  T  U  ■>  N 

4  4*.  (l^lJu)xFfJJCR*X%Y*XPtJP 

1  !) J  H^T(shGTItj*t*XjIofcXtF^«6%F8«6) 

-  1 1  U':!  K 

104  FO  -  Mi.T  t  SnGR  ID*  ,211b,  1P2E16.9,  >H*G, I E/2h*G,l6> 

EM 


Copy  available  to  DTIC  does  not 
Pmnit  iully  legible  reproduction 
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